In thomasgaquiere/Maria: Simulate Different Forms of Logging on a Forest

inventory <- addtreedim(inventorycheckformat(Paracou6_2016), volumeparameters = ForestZoneVolumeParametersTable)

inventory <- treeselection(inventory, objective = 20, scenario ="manual", fuel = "2", diversification = TRUE, specieslax = FALSE,
                           objectivelax = FALSE, topography = DTMParacou, plotslope = PlotSlope, speciescriteria = SpeciesCriteria,
                           advancedloggingparameters = loggingparameters())$inventory

Create a MainTrail and 2 polygons ScndTrail

MainTrail <- sf::st_linestring(matrix(c(286400, 583130,
                                        286400, 583250,
                                        286655, 583250,
                                        286655, 583130,
                                        286400, 583130) # the return
                                      ,ncol=2, byrow=TRUE))

pol1 <- list(matrix(c(286503, 583134,
                      286503, 583240, 
                      286507, 583240, 
                      286507, 583134,
                      286503, 583134) # the return
                    ,ncol=2, byrow=TRUE)) 
pol2 <- list(matrix(c(286650, 583134,
                      286650, 583240, 
                      286654, 583240, 
                      286654, 583134,
                      286650, 583134) # the return
                    ,ncol=2, byrow=TRUE))

PolList = list(pol1,pol2) #list of lists of numeric matrices
ScndTrail = st_multipolygon(PolList)

Succès-échec

# Accessible <- Selected !!!!!!!!!!!! pas oublier pour if (fuel == "0" && directionalfelling != "0")

successfail <- function(
  inventory,
  fuel,
  directionalfelling,
  advancedloggingparameters = loggingparameters()) {

  if (fuel == "0" && directionalfelling == "0"){
    inventory <- inventory %>% 
      mutate(TreeFellingOrientationSuccess = 0)
  }

  if (fuel == "0" && directionalfelling != "0"){
    inventory <- inventory %>% 
      mutate(TreeFellingOrientationSuccess = ifelse(Selected == "1", sample(c(1,0), size = 1, replace = F, prob = c(advancedloggingparameters$TreefallSuccessProportion, 1-advancedloggingparameters$TreefallSuccessProportion)), NA)) # Accessible = linked by 2ndtrails

  }

  if (fuel =="1") {

    inventory <- inventory %>% 
      mutate(TreeFellingOrientationSuccess = ifelse(Selected == "1", sample(c(1,0), size = 1, replace = F, prob = c(advancedloggingparameters$TreefallSuccessProportion, 1-advancedloggingparameters$TreefallSuccessProportion)), NA)) # Selected = not yet linked by 2ndtrails, because 2ndtrails came after

  }

  # Abattre les creux : 
  if (fuel =="2") {

    inventory <- inventory %>% 
      mutate(TreeFellingOrientationSuccess = ifelse(Selected == "1"| ProbedHollow == "1", sample(c(1,0), size = 1, replace = F, prob = c(advancedloggingparameters$TreefallSuccessProportion, 1-advancedloggingparameters$TreefallSuccessProportion)), NA)) # Selected = not yet linked by 2ndtrails, because 2ndtrails came after

  }

  inventory$TreeFellingOrientationSuccess <- as.character(inventory$TreeFellingOrientationSuccess) 

  return(inventory)

}

Rotation function

rotatepolygon <- function(p, angle, fixed) { # angle in the clockwise direction
  p.coords <- sf::st_coordinates(p)[,1:2] # Polygone coordinates extraction
  p.center <- suppressWarnings(sf::st_coordinates(fixed))
  rotate.coords <- function(xy, a, center) {
    co <- cos(-a * pi / 180)
    si <- sin(-a * pi / 180)
    adj <- matrix(rep(center, nrow(xy)), ncol=2, byrow=TRUE) # matrix with fixed point coordinates
    xy <- xy-adj
    cbind(co * xy[,1] - si * xy[,2],si * xy[,1] + co * xy[,2]) + adj
  }
  p.rotate <- rotate.coords(p.coords, a = angle, center = c(p.center[1], 
                                                            p.center[2]))
  Turned <- sf::st_sfc(sf::st_polygon(list(p.rotate))) # create the turned polygon
  return(Turned)
}

Tree creation function

felling1tree <- function(dat,
                    ScndTrail = ScndTrail,
                    MainTrail = MainTrail,
                    fuel, 
                    directionalfelling, 
                    advancedloggingparameters = loggingparameters()){

  Trail <- st_union(MainTrail, ScndTrail) # Our trail will be MainTrail or ScndTrail

  Crown <- dat %>% 
    mutate(xCrown = Xutm, 
           yCrown = Yutm + TrunkHeight + CrownHeight/2,
           exCrown = CrownDiameter/2,
           eyCrown = CrownHeight/2) %>% 
    st_as_sf(coords = c("xCrown", "yCrown")) # ellipse centroid coordinates
  Crown <- st_ellipse(Crown, Crown$exCrown, Crown$eyCrown) # create the ellipse
  Trunk <- with(dat, # and the trunk
                st_polygon(list(matrix(c(Xutm-(DBH/100)/2, Yutm, 
                                         Xutm-(DBH/100)/2, Yutm + TrunkHeight, 
                                         Xutm+(DBH/100)/2, Yutm + TrunkHeight, 
                                         Xutm+(DBH/100)/2, Yutm,
                                         Xutm-(DBH/100)/2, Yutm) # the return
                                       ,ncol=2, byrow=TRUE))))

  # Foot <- st_point(as.numeric(dat[,c("Xutm", "Yutm")])) # Tree coordinates (point)

  # Find the point (TrailPt) on the Trail closest to the location of the tree (Foot)
  Foot <- st_point(c(dat$Xutm,dat$Yutm)) # tree foot point

  NearestPoints <- st_nearest_points(Foot, Trail) # from the Foot of the tree to the Trail (linestring)

  NearestPoint <- st_cast(NearestPoints, "POINT") # to have start (Foot) and end (TrailPt) points
  TrailPt <- NearestPoint[[2]] # the point (TrailPt) on the Trail closest to the location of the tree (Foot)

  theta <- as.numeric(matlib::angle(c(Foot[1] - Foot[1], dat$TreeHeight),
                                    c(TrailPt[1] - Foot[1], TrailPt[2] - Foot[2]), degree = TRUE))

  # Scenarios
  # For a random direction felling
  if(directionalfelling == "0" && (fuel !="1" |fuel !="2")){
    RandomAngle <- sample(c(0:180), size = 1)
    A <- st_difference(st_union(
      rotatepolygon(Trunk, angle = RandomAngle, fixed = Foot), # turned trunk
      rotatepolygon(Crown, angle = RandomAngle, fixed = Foot) # turned crown
    )) 
  }

  # To direct only to avoid damage to future and reserve trees
  # if (directionalfelling == "1"&& (fuel !="1" |fuel !="2")) {
  #   if(dat$TreeFellingOrientationSuccess == "1"){

  # }else{ # else random felling
  # RandomAngle <- sample(c(0:180), size = 1)
  #   A <- st_difference(st_union(
  #   rotatepolygon(Trunk, angle = RandomAngle, fixed = Foot), # turned trunk
  #   rotatepolygon(Crown, angle = RandomAngle, fixed = Foot) # turned crown
  #   ))
  # }
  #   
  # }

  # To direct to avoid damage to future and reserve trees + trail orientation. Winching: Foot before.
  if(directionalfelling == "2"&& (fuel !="1" |fuel !="2")){
    if(dat$TreeFellingOrientationSuccess == "1"){
      A <- st_difference(st_union(
        rotatepolygon(Trunk, angle = 240 + theta, fixed = Foot), # turned trunk
        rotatepolygon(Crown, angle = 240 + theta, fixed = Foot) # turned crown
      ))
      B <- st_difference(st_union(
        rotatepolygon(Trunk, angle = 120 + theta, fixed = Foot), # turned trunk 
        rotatepolygon(Crown, angle = 120 + theta, fixed = Foot) # turned crown
      ))
    }else{ # else random felling
      RandomAngle <- sample(c(0:180), size = 1)

      A <- st_difference(st_union(
        rotatepolygon(Trunk, angle = RandomAngle, fixed = Foot), # turned trunk
        rotatepolygon(Crown, angle = RandomAngle, fixed = Foot) # turned crown
      ))
    }

  }

  # Fuel wood exploitation in the crowns
  if(fuel =="1" |fuel =="2"){

    TrailDist <- st_distance(Foot, TrailPt)

    # ADD SLOPE CRITERIA !!!

    if(dat$TreeFellingOrientationSuccess == "1"){

      if(TrailDist <= advancedloggingparameters$GrappleLength){ # <= 6m (= grapple length) -> winching by grapple
        A <- st_difference(st_union(
          rotatepolygon(Trunk, angle = theta + 60 , fixed = Foot), # turned trunk
          rotatepolygon(Crown, angle = theta + 60 , fixed = Foot) # turned crown
        ))
        B <- st_difference(st_union(
          rotatepolygon(Trunk, angle = 300 + theta, fixed = Foot), # turned trunk
          rotatepolygon(Crown, angle = 300 + theta, fixed = Foot) # turned crown
        ))
      } else { # > 6m -> winching by cable
        A <- st_difference(st_union(
          rotatepolygon(Trunk, angle = 240 + theta, fixed = Foot), # turned trunk
          rotatepolygon(Crown, angle = 240 + theta, fixed = Foot) # turned crown
        ))
        B <- st_difference(st_union(
          rotatepolygon(Trunk, angle = 120 + theta, fixed = Foot), # turned trunk 
          rotatepolygon(Crown, angle = 120 + theta, fixed = Foot) # turned crown
        ))
      }
    }else{ # else random felling
      RandomAngle <- sample(c(0:180), size = 1)
      A <- st_difference(st_union(
        rotatepolygon(Trunk, angle = RandomAngle, fixed = Foot), # turned trunk
        rotatepolygon(Crown, angle = RandomAngle, fixed = Foot) # turned crown
      ))
    }

  }

  FellingOuputs <- list(Foot = Foot,
                        NearestPoints = NearestPoints,
                        TrailPt = TrailPt,
                        A = A)

  return(FellingOuputs)

}

Application

# Compute fellingtree succes and fails
inventory <- successfail(
  inventory,
  fuel = "2",
  directionalfelling = "2",
  advancedloggingparameters = loggingparameters())

# Only selected trees
SelectedTrees <- inventory %>% 
  filter(Selected == "1") %>%
  dplyr::select(idTree,DBH,TrunkHeight,TreeHeight,CrownHeight,CrownDiameter,Selected, Xutm, Yutm, TreeFellingOrientationSuccess)

# To extract de the shortest ways
NearestPoints <- SelectedTrees %>%
  group_by(idTree) %>% 
  do(geometry = felling1tree(., ScndTrail = ScndTrail, MainTrail = MainTrail, fuel = "2", directionalfelling ="2", advancedloggingparameters = loggingparameters())$NearestPoints) %>% 
  mutate(geometry = as(geometry, "sfc_LINESTRING")) %>% # list -> sfc_MULTIPOLYGON
  st_as_sf()

# Application in a dataframe
FellingTrees <- SelectedTrees %>% 
  # mutate(angle = list(c(0, 30))) %>% # for each row, a list of 2 values -> a dataframe in a dataframe
  # unnest(angle) %>% # display the small dataframe in the big dataframe in one
  group_by(idTree) %>% # for each tree and each given angle 
  do(geometry = felling1tree(., ScndTrail = ScndTrail, MainTrail = MainTrail, fuel = "2", directionalfelling ="2", advancedloggingparameters = loggingparameters())$A) %>% # inform geometry. # Filling a column from a function whose input is a table


  # or version more Sf (to test)
  # pol <- st_as_sf(SelectedTrees, 
  #                 as(geometry = st_as_sfc(
  #                   lapply(split(SelectedTrees, SelectedTrees$idTree), felling1tree(30))
  #                 )
  #                 , "sfc_MULTIPOLYGON")
  # )



library(ggplot2)
ggplot() +
  geom_sf(data = st_union(Trunk, Crown), colour = "red") +
  geom_sf(data = Foot) +
  geom_sf(data = Trail) +
  geom_sf(data = TrailPt) +
  geom_sf(data = NearestPoints) +
  geom_sf(data = A)+
  geom_sf(data = B)

FellingTrees <- inventory %>% 
  group_by(idTree) %>% # for each tree and each given angle 
  do(geometry =  # inform geometry. # Filling a column from a function whose input is a table
       ifelse(!is.na(.$TreeFellingOrientationSuccess),
              felling1tree(.,
                      ScndTrail = ScndTrail,
                      MainTrail = MainTrail,
                      fuel = "2", directionalfelling ="2",
                      advancedloggingparameters = loggingparameters())$A,
              NULL)) %>%
  mutate(geometry = as(geometry, "sfc_MULTIPOLYGON")) %>% # list -> sfc_MULTIPOLYGON
  st_as_sf()

felttress <- NA
inventory %>% 
  filter(!is.na(TreeFellingOrientationSuccess)) %>% 
  group_by(idTree) %>% # for each tree and each given angle 
  do(geometry =  # inform geometry. # Filling a column from a function whose input is a table
       felling1tree(.,
               ScndTrail = ScndTrail,
               MainTrail = MainTrail,
               fuel = "2", directionalfelling = "2",
               advancedloggingparameters = loggingparameters())$A) %>%
  mutate(geometry = as(geometry, "sfc_MULTIPOLYGON")) %>% # list -> sfc_MULTIPOLYGON
  st_as_sf() %>% 
  mutate(geometry = st_as_text(geometry))


left_join(inventory, poly) %>% 
  ggplot() +
  geom_sf()

st_as_text(poly$geometry)


felttrees <- inventory %>% 
  filter(!is.na(TreeFellingOrientationSuccess)) %>% 
  felling1tree(.,
          ScndTrail = ScndTrail,
          MainTrail = MainTrail,
          fuel = "2", directionalfelling = "2",
          advancedloggingparameters = loggingparameters())$A %>% 
  st_as_text()  # as text to easy join with a non spacial table
# unnest(treepolygon) # display the small dataframe in the big dataframe in one

inventory <- left_join(inventory, felttrees) # join spatial filtered inventory and non spatial complete inventory

get_geometry <- function(inventory){ # function from text (wkt) to sf
  st_as_text() %>% 
    unnest(treepolygon)
}

inventory2 <- left_join(inventory, felttrees)

get_geometry <- function(inventory){
  inventory %>% 
    filter(!is.na(treepolygon)) %>%
    st_as_sf(wkt = "treepolygon")
}

ggplot() +
  geom_sf(data = st_as_sf(inventory, coords = c("Xutm", "Yutm"))) + # all the individuals of the plot
  geom_sf(data = get_geometry(inventory), fill = "red") # trees polygons

st_intersection( # trees under the fallen trees
  get_geometry(inventory),
  geom_sf(data = st_as_sf(inventory, coords = c("Xutm", "Yutm"))) +
    geom_sf(data = get_geometry(inventory2), fill = "red"))

st_intersection(
  get_geometry(inventory2),
  st_as_sf(inventory, coords = c("Xutm", "Yutm"))
) %>% 
  ggplot() +
  geom_sf()

inventory <- successfail(
  inventory,
  fuel = "2", directionalfelling = "2",
  advancedloggingparameters = loggingparameters())

# Treefelling
felttrees <- inventory %>%
  filter(!is.na(TreeFellingOrientationSuccess)) %>%
  group_by(idTree) %>% # for each tree
  do(TreePolygon = # inform geometry. # Filling a column from a function whose input is a table
       felling1tree(.,
               fuel = "2", directionalfelling = "2",
               MainTrail = MainTrail, ScndTrail = ScndTrail,
               advancedloggingparameters = loggingparameters())$A %>%
       st_as_text()) %>% # as text to easy join with a non spacial table
  tidyr::unnest(TreePolygon) # here to pass from list to character

inventory <- left_join(inventory, felttrees, by = "idTree") # join spatial filtered inventory and non spatial complete inventory